rm(list = ls())
library(ggplot2)
library(lme4)
library(BayesFactor)
library(plyr)
setwd('~/Dropbox/WoMAAC/Experiments/UK CH/Phase 2/xInt_1_2_3/')
xInt <- read.csv('P2_YA_xInt.csv')

logit2odds <- function(logit){
  tmp <- logit
  if(logit<1){
    tmp <- tmp*-1}
  tmp <- tmp*2
  return(exp(tmp))
}

SAVE = F

xInt$ID <- paste(xInt$participant, xInt$experiment, sep = '_')
xInt$ID <- as.factor(xInt$ID)

xInt$memDomain[xInt$experiment %in% c('MemVerb_ProcVerb', 'MemVerb_ProcVis')] <- 'Verbal'
xInt$memDomain[xInt$experiment %in% c('MemVis_ProcVerb', 'MemVis_ProcVis')] <- 'Visual'
xInt$procDomain[xInt$experiment %in% c('MemVerb_ProcVerb', 'MemVis_ProcVerb')] <- 'Verbal'
xInt$procDomain[xInt$experiment %in% c('MemVerb_ProcVis', 'MemVis_ProcVis')] <- 'Visual'
xInt$memDomain <- as.factor(xInt$memDomain)
xInt$procDomain <- as.factor(xInt$procDomain)

xInt <- xInt[xInt$practice == 0,]
mainData_mem <- xInt[xInt$condition %in% c('ST_Mem', 'DT'),]
mainData_mem <- droplevels(mainData_mem)
mainData_mem$condition <- relevel(mainData_mem$condition, ref = 'ST_Mem')

mainData_proc <- xInt[xInt$condition %in% c('ST_Proc', 'DT'),]
mainData_proc <- droplevels(mainData_proc)
mainData_proc$condition <- relevel(mainData_proc$condition, ref = 'ST_Proc')

### Spans

spans <- ddply(xInt, c('ID', 'experiment', 'memDomain', 'procDomain'), summarize,
              memSpan = mean(memSpan, na.rm = T), procSpan = mean(procSpan, na.rm = T))

spans <- within(spans, {
  ID <- as.factor(ID)
  experiment <- as.factor(experiment)
  memDomain <- as.factor(memDomain)
  procDomain <- as.factor(procDomain)
})

with(spans, table(experiment))

anovaBF_memSpan <- anovaBF(memSpan ~ experiment*memDomain, data = spans)
anovaBF_memSpan[2]/anovaBF_memSpan[3] # BF = 8.36 in favour of memDomain only

anovaBF_procSpan <- anovaBF(procSpan ~ experiment*procDomain, data = spans)
anovaBF_procSpan[1]/anovaBF_procSpan[3] # BF = noninformative

contrasts(mainData_mem$condition) <- c(-1,1)
contrasts(mainData_mem$memDomain) <- c(-1,1)
contrasts(mainData_mem$procDomain) <- c(-1,1)

# Memory analysis
m_full_mem <- glmer(cbind(MemNcorr, memSpan - MemNcorr)
                             ~1 + condition*memDomain*procDomain
                             +(1|ID),
                             data = mainData_mem,
                             family = binomial(link = logit)
)

# test removal of three way interaction
m_1_mem <- update(m_full_mem, .~. - condition:memDomain:procDomain)
(bic_mem_fullvs1 <- BIC(m_full_mem) - BIC(m_1_mem)) # 4.57, BF = 9.84 for removal

# test removal of two way interactions
m_2_mem <- update(m_1_mem, .~. - memDomain:procDomain)
m_3_mem <- update(m_1_mem, .~. - condition:procDomain)
m_4_mem <- update(m_1_mem, .~. - condition:memDomain)

# memDomain:procDomain
(bic_mem_1vs2 <-  BIC(m_1_mem) - BIC(m_2_mem)) # 8.35, BF = 65.18 for removal
#condition:procDomain
(bic_mem_1vs3 <- BIC(m_1_mem) - BIC(m_3_mem)) # 8.40, BF = 66.57 for removal
#condition:memDomain
(bic_mem_1vs4 <- BIC(m_1_mem) - BIC(m_4_mem)) # -106.80, BF = massive for retention

# remove condition:memDomain interaction
m_5_mem <- update(m_1_mem, .~. - memDomain:procDomain - condition:procDomain)

# test removal of main effect of procDomain
m_6_mem <- update(m_5_mem, .~. - procDomain)
(bic_mem_5vs6 <- BIC(m_5_mem) - BIC(m_6_mem)) # 8.40, BF = 66.56 for removal

# Winning model
m_final_mem <- m_6_mem
summary(m_final_mem)

# confidence intervals
ci <- confint(m_final_mem)

ci_table_mem <- data.frame(cbind(fixef(m_final_mem), ci[2:5,]))

names(ci_table_mem) <- c('coeff', 'lower', 'upper')

ci_odds_mem <- apply(ci_table_mem, c(1,2), logit2odds)


# Processing analysis
contrasts(mainData_proc$condition) <- c(-1,1)
contrasts(mainData_proc$procDomain) <- c(-1,1)
contrasts(mainData_proc$memDomain) <- c(-1,1)

m_full_proc <- glmer(cbind(procNcorr, procSpan - procNcorr)
                       ~1 + condition*memDomain*procDomain
                       +(1|ID),
                       data = mainData_proc,
                       family = binomial(link = logit)
)

# test removal of three way interaction
m_1_proc <- update(m_full_proc, .~. - condition:memDomain:procDomain)
(bic_proc_fullvs1 <- BIC(m_full_proc) - BIC(m_1_proc)) # 8.10, BF = 57.40 to remove

# test removal of two way interactions
m_2_proc <- update(m_1_proc, .~. - memDomain:procDomain)
m_3_proc <- update(m_1_proc, .~. - condition:procDomain)
m_4_proc <- update(m_1_proc, .~. -condition:memDomain)

# memDomain:procDomain
(bic_proc_1vs2 <- BIC(m_1_proc) - BIC(m_2_proc)) # 8.18, BF = 59.60 to remove
# condition:procDomain
(bic_proc_1vs3 <- BIC(m_1_proc) - BIC(m_3_proc)) # 5.74, BF = 17.59 to remove
# condition:memDomain
(bic_proc_1vs4 <- BIC(m_1_proc) - BIC(m_4_proc)) # 8.40, BF = 66.68 to remove

# remove all three interactions
m_5_proc <- update(m_1_proc, .~. - condition:memDomain - condition:procDomain - memDomain:procDomain)

# test removal of main effects
m_6_proc <- update(m_5_proc, .~. - procDomain)
m_7_proc <- update(m_5_proc, .~. - memDomain)
m_8_proc <- update(m_5_proc, .~. - condition)

# procDomain
(bic_proc_5vs6 <- BIC(m_5_proc) - BIC(m_6_proc)) # 6.22, BF = 22.37 to remove
# memDomain
(bic_proc_5vs7 <- BIC(m_5_proc) - BIC(m_7_proc)) # 6.85, BF = 30.72 to remove
# condition
(bic_proc_5vs8 <- BIC(m_5_proc) - BIC(m_8_proc)) # -90.86, BF = massive to retain

# winning model
m_final_proc <- update(m_5_proc, .~. - procDomain - memDomain)
(bic_proc_5vsfinal <- BIC(m_5_proc) - BIC(m_final_proc))
summary(m_final_proc)

# confidence intervals
ci <- confint(m_final_proc)

ci_table_proc <- data.frame(cbind(fixef(m_final_proc), ci[2:3,]))

names(ci_table_proc) <- c('coeff', 'lower', 'upper')

ci_odds_proc <- apply(ci_table_proc, c(1,2), logit2odds)


if (SAVE){
  dir.create('Rdata_files/', showWarnings = F)
  
  save(list = c(ls(pattern = 'm_'), ls(pattern = 'bic_'), ls(pattern = 'anovaBF'), ls(pattern = 'ci_') , 'mainData_mem', 'mainData_proc', 'spans'), file = 'Rdata_files/P2_YA_main_xInt.RData')
  
}
